Brownian motion and Diffusion
Contents
Brownian motion and Diffusion#
import numpy as np
import scipy
from numpy.random import normal, choice, uniform
import ipywidgets as widgets
import matplotlib.pyplot as plt
plt.style.use('fast')
%matplotlib inline
Mean square displacement (MSD) of random walker#
After time t number of steps (time) how far has random walker moved from the origin?
We quantify this by computing means square displacement. Where mean is computed over N number of simulated (observed) trajectories.
Einstein developed a theory of diffusion based on random walk ideas and obtained a key equation relating mean square displacement to time \((t = n \Delta t)\)
MSD of a 1D random walk#
def rw_1d(T, N):
'''
L: trajectory length
N: Number of trajecotries
returns np.array with shape (T, N)
'''
# Create random walks
r = choice([-1,1], size=(T, N))
#Accumulate position
rw = r.cumsum(axis=0)
#Set initial position
rw[0,:]=0
return rw
T, N = 2000, 1000
rw = rw_1d(T, N)
t = np.arange(T)
R2 = (rw[:, :] - rw[0, :])**2 # Notice we subtract initial time
msd = np.mean(R2, axis=1) # Notice we average over N
plt.loglog(t, np.sqrt(msd), lw=3)
plt.loglog(t, np.sqrt(t), '--')
plt.title('Compute mean square deviation of 1D random walker',fontsize=15)
plt.xlabel('Number of steps, n',fontsize=15)
plt.ylabel(r'$MSD(n)$',fontsize=15);
MSD of a 2D random walk#
def rw_2d(T, N):
'''2d random walk function:
T: trajectory length
N: Number of trajecotry
returns np.array with shape (T, N)
'''
verteces = np.array([(1, 0),
(0, 1),
(-1, 0),
(0, -1)])
rw = verteces[choice([0,1,2,3], size=(T, N))]
rw[0, :, :] = 0
return rw.cumsum(axis=0)
traj = rw_2d(T=10000, N=1000)
#Simulate 2D random walk
T, N = 10000, 1000
traj = rw_2d(T, N)
#Compute RSD
dx = (traj[:, :, 0]- traj[0, :, 0])
dy = (traj[:, :, 1]- traj[0, :, 1])
R2 = np.mean(dx**2 + dy**2, axis = 1) # notice how we averaging over N
fig, ax = plt.subplots(nrows=2, figsize=(10,10))
t = np.arange(T) # time axis
ax[1].loglog(t, np.sqrt(R2), lw=3, alpha=0.5);
ax[1].loglog(t, np.sqrt(t), '--');
ax[0].set_title('2D random walker',fontsize=15)
ax[0].plot(traj[:3000, :5, 0], traj[:3000, :5, 1]);
ax[0].set_xlabel('X')
ax[1].set_xlabel('Y')
ax[1].set_xlabel('Number of steps, n',fontsize=15)
ax[1].set_ylabel(r'$MSD(n)$',fontsize=15);
Brownian motion#
The Brownian motion describes the movement of a particle suspended in a fluid resulting from random collisions with the quick molecules in the fluid (diffusion).
A small particle undergoes large number of molecular collisions when going from one step to another or after \(dt\) time. As a result each displacement over time \(dt\) can be viewed as a sum of ranomd collisions which can be approximated by a normal distribution via Central Limit Theorem.
Thus to simulate brownian motion we draw ranom displacements from normally distribution.
We assume we have started at position \(\mu=0\) and our variance is given by \(\sigma^2=2Dt\) Where D is the diffusion coefficnets which is related to parameters of discree random walk as shown in the lecture.
In the last step we re-wrote brownian motion equation in a convenient way by shifting normally distributed radnom variable by \(\mu\) and scaling it by \(\sigma\)
def brown(T, N, dt=1, D=1):
"""
Creates 3D brownian path given:
time T
N=1 trajecotires
dt=1 timestep
D=1 diffusion coeff
returns np.array with shape (N, T, 3)
"""
n = int(T/dt) # how many points to sample
dR = np.sqrt(2*D*dt) * np.random.randn(N, n, 3) # 3D position of brownian particle
R = np.cumsum(dR, axis=1) # accumulated 3D position of brownian particle
return R
Below we proceed to simulate continuous random walk in 1D-3D
We will plot trajecotires and distributions of brownian particle using interactive plotting via ipywidgets and holoviews/plotly interface.
R=brown(T=3000, N=1000)
print(R.shape)
(1000, 3000, 3)
@widgets.interact(t=(10,3000-1))
def brownian_plot(t=10):
fig, ax = plt.subplots(ncols=2)
ax[0].plot(R[:5, :t, 0].T, R[:5, :t, 1].T);
ax[1].hist(R[:, 10, 0], density=True, color='red')
ax[1].hist(R[:, t, 0], density=True)
ax[1].set_ylim([0,0.1])
ax[0].set_ylim([-200, 200])
ax[0].set_xlim([-200, 200])
fig.tight_layout()
import holoviews as hv
hv.extension('plotly')
plots = []
for i in range(10):
plot = hv.Path3D(R[i,:,:], label='3D random walk').opts(width=600, height=600, line_width=5)
plots.append(plot)
hv.Overlay(plots)
rw_curve = [hv.Curve((R[i,:,0], R[i,:,1])) for i in range(10)]
xdist = hv.Distribution(R[:,10,0], ['X'], ['P(X)'])
ydist = hv.Distribution(R[:,10,1], ['Y'], ['P(Y)'])
hv.Overlay(rw_curve) << ydist << xdist
Diffusion Equation#
The movement of individual random walkers \(\leftrightarrow\) density of walkers \(\rho(\vec{r},t)\)
Diffusion equation:
Formulated empirically as Fick’s laws
This is a 2nd order PDE! Unlike equations of motion diff eq shows irreersibile behaviour
This one exactly solvable. But in general reaction-diffusion PDEs difficult to solve analytically.
Can solve numerically by writing derivatives as finite differences!
Can also simulate via random walk!
Diffusion coefficient \(D\), Units \([L^2]/[T]\)
Important special case solution (here written in 1d):
where \(\sigma(t)=\sqrt{2{D}t}\)
density remains Gaussian
Gaussian becomes wider with time
check that this is indeed a solution by plugging into the diffusion equation!
def sigma(t, D = 1):
return np.sqrt(2*D*t)
def gaussian(x, t):
return 1/np.sqrt(2*np.pi*sigma(t)**2) * np.exp(-x**2/(2*sigma(t)**2)) #
@widgets.interact(t=(1,100,1))
def diffusion(t=0.001):
R = brown(T=101, N=1000)
x = np.linspace(-20, 20, 100)
plt.plot(x, gaussian(x, 1), '--', color='orange', label='t=0')
plt.plot(x, gaussian(x, t), lw=3, color='green', label=f't={t}')
plt.hist(R[:,t,0], density=True, alpha=0.6, label='simulation hist')
plt.legend()
plt.ylabel('$p(x)$')
plt.xlabel('$x$')
plt.xlim([-25, 25])
References#
The mighty little books
More in depth
On the applied side